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Abstract 

We propose a proximal approach to deal with convex optimization problems involving nonlinear 
constraints. A large family of such constraints, proven to be effective in the solution of inverse problems, 
can be expressed as the lower level set of a sum of convex functions evaluated over different, but 
possibly overlapping, blocks of the signal. For this class of constraints, the associated projection operator 
generally does not have a closed form. We circumvent this difficulty by splitting the lower level set into 
as many epigraphs as functions involved in the sum. A closed half-space constraint is also enforced, in 
order to limit the sum of the introduced epigraphical variables to the upper bound of the original lower 
level set. 

In this paper, we focus on a family of constraints involving linear transforms of £i tP balls. Our main 
theoretical contribution is to provide closed form expressions of the epigraphical projections associated 
with the Euclidean norm (p = 2) and the sup norm (p — +00). The proposed approach is validated in the 
context of image restoration with missing samples, by making use of TV-like constraints. Experiments 
show that our method leads to significant improvements in term of convergence speed over existing 
algorithms for solving similar constrained problems. 



1 Introduction 

As an offspring of the wide interest in frame representations and sparsity promoting techniques for data 
recovery, there has been an emergence of proximal methods to efficiently solve convex optimization problems 
arising in inverse problems. Examples of applications of these methods can be found in f-MRI reconstruction 
[TJ [2] , satellite image restoration [3J H] , microscopy image deconvolution [SJ |5] , computed tomography [7] , 
Positron Emission Tomography [5J [5], texture-geometry decomposition [IUJ [TTJ [T2] , machine learning [T31[T1] . 
stereo vision [TS], and audio processing [T^ITT]. Proximal algorithms have gained much popularity in solving 
large-size optimization problems involving non-diffcrcntiable (or even non finite) functions. One of the main 
advantages of these methods is that they are amenable to parallel implementations. For a survey on proximal 
algorithms and their applications, the reader is referred to |T8j [13] . Note also that some of these methods 
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are closely related to augmented Lagrangian approaches (THl [2S1 HI] ■ Even if proximal algorithms and the 
associated convergence properties have been deeply investigated [22J [23l [24J [25] , some questions persist in 
their use for solving inverse problems. A hrst question is: how can we set the parameters serving to enforce 
the regularity of the solution in an automatic way? Another question is related to the selection of the most 
appropriate algorithm within the class of proximal algorithms for a given application. This also raises the 
question of the computation of the proximity operators associated with the different functions involved in 
the criterion. Various strategies were proposed in order to address the first question [26] [27] [28l [29] [30] , but 
the computational cost of these methods is often high, especially when several regularization parameters 
have to be set. Alternatively, it has been recognized for a long time that incorporating constraints directly 
on the solutions [32 [Ml 1331 1311 135] instead of considering regularized functions may often facilitate the 
choice of the involved parameters. Indeed, in a constrained formulation, the constraint bounds are usually 
related to some physical properties of the target solution or some knowledge of the degradation process, 
e.g. the noise statistical properties. Note also that there exist some conceptual Lagrangian equivalences 
between regularized solutions to inverse problems and constrained ones, although some caution should be 
taken when the regularization functions are nonsmooth (see [36] where the case of a single regularization 
parameter is investigated). 

In this context, the objective of this paper is to provide an answer to the second question, that is to 
propose efficient parallel algorithms for solving the following constrained convex optimization problem: 

Problem 1.1. 

R 

minimize g r (T r x) s. t. < ■ (1) 

r=l 

where 

(i) for every re {1, . . . , R}, T r is a bounded linear operator from a real Hilbert space T-L to M. Nr ; 

(ii) for every r € {1, . . . , R}, g r is a proper lower-semicontinuous convex function from M. Nr to ]— oo, +oo]; 
(hi) for every s G {1, . . . , S}, H s is a bounded linear operator from T-L to M A/s ; 

(iv) for every s € {1, . . . , S}, C s is a nonempty closed convex subset of R Ms . 

A large number of inverse problems can be formulated under the form of Problem 11.11 A classical 
application example is image recovery from blurred and noisy observations. Let z denote the vector of 
observed data. When the noise is assumed to be zero-mean additive white Gaussian, the restored data can 
be obtained by solving Problem 11.11 where R = 2, S = 0, g\ = A|| • ||i with A € ]0, +oo[, T\ is an analysis 
frame operator [37] allowing us to sparsify signal x, g 2 = || • — z\\ 2 , and T 2 denotes the matrix associated 
with the degradation blur [3H1 H3] • On the other hand, a constrained formulation of the same restoration 
problem leads to Problem 11.11 with R = 1, S = 1, g\ = || • ||i, T\ is the same frame operator as previously, 
and 

(VieH) HxxeCx <S> \\H lX -z\\ 2 < ?7i. 

Here, H± denotes the matrix associated with the degradation blur, and r\\ is a positive constant which is 
typically chosen proportional to the noise variance. This specific constraint formulation has been considered 
in [35]. 

The present work aims at designing efficient methods in order to deal with Problem 11.11 when the convex 
constraints are expressed as follows: for every s € {1, . . . , S}, 

(Vz e H) H s x eC s & h s (H s x) < r? s , (2) 
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where rj s £E K and h s is a proper lower-semicontinuous function from R Ms to ]— oo, +00]. 

The paper is organized as follows. In Section [21 we motivate the choice of proximal tools and recall 
some of their theoretical properties. In addition, closed form expressions for specific proximity operators are 
derived. Then, in order to deal with a constraint expressed under the general form ([2]) , a splitting approach 
involving an epigraphical projection is proposed. This epigraphical projection technique is described in 
detail in Section [3] Experiments in the context of image reconstruction are presented in Section 2) Finally, 
some conclusions are drawn in Section 

Notation: Let H be a real Hilbert space endowed with the norm || • || and the scalar product (• | ■). T a {H) 
denotes the set of proper lower-semicontinuous convex functions from % to ]— 00, +00]. Recall that a function 
ip: T~L — > ]— 00, +00] is proper if its domain domip = {y 6 T-L | <p(y) < +00} is nonempty. The epigraph of 
ip G To('H) is the nonempty closed convex subset of^xl defined as epi^ = {(t/, £) e % x R | <p(y) < 
and the lower level set of ip at height £ € M is the nonempty closed convex subset of H defined as lev<^ (p = 
{y £ W. I <p(y) < (}. A subgradient of ip at y G W is an element of its subdifferential defined as d<p(y) = 
{t £ % J (Vu € H) >p(u) > tp(y) + (t I u — y)}Q The indicator function l c € T (H) of a nonempty closed 
convex subset C of % is given by 

(VyeH) , c (y) = \°; \^ C : (3) 

l+oo, otherwise. 
The relative interior of a subset C of H is denoted by ri C. 



2 Proximal tools 

2.1 Prom gradient descent to proximal algorithms 

The first methods for finding a solution to an inverse problem were restricted to the use of a differentiable cost 
function [40], i.e. Problem lf.il where S = and, for every r £ {1, . . . , R}, g r denotes a differentiable function. 
In this context, gradient-based algorithms, e.g. nonlinear conjugate gradient or quasi-Newton methods, are 
popular (see |3T] and references therein for recent developments concerning these approaches). However, in 
order to model additional properties, sparsity promoting penalizations or hard constraints (S > 1) may be 
introduced and the diffentiability property is not satisfied anymore. One way to circumvent this difficulty is 
to resort to smart approximations in order to smooth the involved non-differentiable functions [421 431 144) . 
If one wants to address the original nonsmooth problem without introducing approximation errors, one may 
apply some specific algorithms e.g. Gauss-Seidel or Uzawa methods, the convergence of which is guaranteed 
under restrictive assumptions [45]. Interior point methods [46] can also be employed for small to medium 
size optimization problems. 

On the other hand, in order to solve convex feasibility problems, i.e. to find a vector belonging to the 
intersection of convex sets (Problem 11.11 with R — 0), iterative projection methods were developed. The 
projection onto convex sets algorithm (POCS) is one of the most popular approach to solve data recovery 
problems [371 HH1 GUI HH]- A drawback of POCS is that it is not well-suited for parallel implementations. 
The Parallel Projection Method (PPM) and Method of Parallel Projections (MOPP) are variants of POCS 
making use of parallel projections. Moreover, these algorithms were designed to efficiently solve inconsistent 
feasibility problems (when the intersection of the convex set is empty). Thorough comparisons between 
projection methods have been performed in |50[ 151] . 

1 When ip is Gatcaux-differcntiable at y, d<p(y) = {V(/?(j/)} where \7ip(y) is the gradient of ip at y. 
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Computing the projection Pc onto a nonempty closed convex subset C of a real Hilbert space TL requires 
to solve a constrained quadratic minimization problem: 

(Vy e 'H) Pc{y) = argmin||u -y\\. (4) 

The distance to C of every point y G TL is then given by dc(y) = \\y — Pc(y)\\- However, it turns out 
that a closed form expression of the solution to (U) is available in a limited number of instances. One such 
well-known example is the projection onto a closed half-space: 

Proposition 2.1. [HJ Lei 77 G M, fe£ f G H\{0}, and let C = {u € H\ (u\t) < n}. The projection ofy G % 
onto C is expressed as 

[y, if (y\t)<v, 

C ^ I y + - - t, otherwise. ^ 

k \\t\r 



When an expression of the direct projection is not available, the convex set C can be approximated by 
a half-space, which leads to the concept of subgradient projection. 

Proposition 2.2. [44] Let 77 G K and let ip £ To('H). Suppose that C = lev< n ip ^ 0. Lei y E W and let 

t £ dip(y) be a subgradient of ip at y. Then, C is a subset of the half-space 

Cy = {ueH I (u-y\t) <r,-<p(v)}. (6) 
Ifty^O, the projection ofy onto is a subgradient projection of y G "H onto C. 



The subgradient projection plays a key role in Polyak algorithm that alternates a subgradient projection 
and an exact projection [52] (see also [53] for alternative projected subgradient approaches). An efficient 
block iterative surrogate splitting method was proposed in |54) in order to solve Problem ll.il when, for every 
rE {!,..., R], g r = || • -z r \\ 2 where z r G M. Nr . A main limitation of this method is that the global objective 
function must be strictly convex. For recent works about subgradient projection methods, the readers may 
refer to |55l [56] ■ 



A way to overcome this difficulty consists of considering proximal approaches. The key tool in these 
methods is the proximity operator [57 of a function p € To('H), defined as 

{VyETt) prox (y) = argmini||u - y\\ 2 + ip{u). (7) 

The proximity operator generalizes the notion of projection onto a nonempty closed convex subset C of Ti 
in the sense that prox tc = Pc- For every y G Ti. p = prox ¥ ,(y) is uniquely defined through the inclusion 

y-pedp(p). (8) 

Proximity operators enjoy many additional interesting properties. Some of them are recalled next. 

Example 2.3. [23] Let t G ]0,+oo[, (3 G [l,+oo[, and set 

ip: E -> ]-oo,+oo] : ^ T \^f . (9) 

Then, for every £ G M, prox^^) is given by 

'sign(f) max{|£| - r, 0}, 



1+2-r' 



lsign(0* 67 > 





1, 


ifP = 


4 
3 


ifP = 


3 
2 


if P = 


2, 


ifP = 


3. 



i 
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where e = ^ 2 + 256r 3 /729 and sign is the signum function. 



It can be noticed that the proximity operator associated with j5 = 1 reduces to a soft thresholding [58] . 

The class of proximal methods includes primal algorithms [5DJ [221 (2H EPJ EH 1E21 ESS EEl HI] and primal- 
dual algorithms [Ml EUS [DDJ EZJ [DH1 (Ml ED]. Primal algorithms generally require to compute inverses of 
some linear operators (typically, X2r=i T*T r + X)f=i ^s^s), while primal-dual ones only require to com- 
pute (T r )!< r <fl and (i? s )i< s <5, and their adjoints (T*)i< r <R and (H*)i< s <s- Consequently, primal-dual 
methods are often easier to implement than primal ones, but their convergence may be slower [711 172) . 



2.2 Proximity operators: new closed forms 

The projection onto a convex set such as defined in IJ5J often is non trivial. In Section [3[ we will show that 
this problem can be solved by resorting to a set of epigraphical projections which are easier to compute. 
The key point is that epigraphical projections are closely related to proximity operators. In the following, 
we provide some results that will be useful for the calculation of the projection onto the epigraph of a convex 
function (all the proofs are provided in the appendix). 

Proposition 2.4. Let % be a real Hilbert space and letHxM. be equipped with the standard product space 
norm. Let if be a function in To('H) such that doimp is open. The projector P epiv onto the epigraph of ip 
is given by: 

(V(yX)eHxR) P Bpiv to,C) = (p,9) (11) 

where 

ip =proxi (roax{v ,_ f)0})2 (y), 
[6 = max{^(p),C}. 

Note that alternative characterizations of the epigraphical projection can be found in J73J Proposi- 
tions 9.17, 28.28]. 

Let tp be a function in Fo(H) such that dom^s is open. From Proposition 12 .41 we see that 
proxi( max j <(3 _^ i\2 with £ G R plays a prominent role in the calculation of the projection onto epiip. We 
now provide examples of functions if for which this proximity operator admits a simple form. 

Proposition 2.5. Let (3 £ [1, +oo[, t G ]0, +oo[. Assume that 

(Vy G R) <p(y) = r\yf. (13) 

If £ g] — oo, 0], then for every y G R 

r sign(y) if/3 = l, .... 

Proxi (max{ ^ co})2 (y) - <^ 1 + r 2 J ' (14) 

[sign(y)xo, if P > 1, 

where xo is the unique solution on [0, +oo[ o/ i/ie equation 

Pr 2 x 2 ^ 1 -PTCx^ 1 +X=\y\- (15) 

If C £ ]0; i^en, /or every y G R, 

prox l(maxW , 0})2(y) = { y : */M<(cAr, (16) 

[sign(j/)x (c/r )i/^, otherwise, 



where Xu/tY/p * s ^ e unique solution on [(■|) 1 ^' 3 , +oo[ o/ JT5 
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Note that, when /3 is a rational number, (|15[) is equivalent to a polynomial equation for which either 
closed form solutions are known or standard numerical solutions exist. 

Proposition 2.6. Let % be a real Hilbert space and let C be a nonempty convex subset ofH. Let j3 G [1, +00 [, 

r G ]0, +00 [, and (el. Then, for every y G H, 

( x _ jy, ify e c, 

P««i(max { r4-c,o})»W " i ay + (1 _ a ) Pc ^ otherwise, { 7 > 



P rOX i(max{r|.|/3-C,0})2 (<fc(2/)) 

re a = 

1231 



where a = 2 max T ^ - ' and £/ie expression o/proxi ( max { T |.|/s_f o}) 2 * s provided by Proposition 



The following result is a consequence of the previous one: 
Corollary 2.7. Let % be a real Hilbert space. Assume that 

(VyeH) <p(y)=T\\y-z\\ (18) 
where z g H, ( £ I and r G ]0, +oo[. Then, for every i/gK, 

!z, if \\y-z\\ < ~ T C, 

y, if\\y-A\<^ ( 19 ) 

z + a(y — z), otherwise, 

where a = ^{l + j^ 

Another result which will be used is the following: 
Proposition 2.8. Let M G N* and let <p be defined as 

1 M 

{\fy G R) V(») = 5 E (max{r( m )(^-y),0}) 2 (20) 

rn — 1 

i/Aere (r«, . . . , T ( M )) T G R M , and (i/W, . . . , z/ M )) T G M M is sucft i/mi < ... < v {M) . Set = -oo 
and j/( M+1 ) = +oo. Then, ip G To(K) 7 and for every y G M., 

where, for every m G {1, . . . , M }, t!"^ = min(T ( ™- ) , 0) and ri_ = max(r^ m ^ , 0), and m is the unique integer 
in {1, . . . , M + 1} smc/i £/ia£ 



i+E:^(^ m) ) 2 +E^(4 m) ) 2 



(wii/i ifte convention Em=i ' = Em=Af+i ' = 



G 



3 Projection computation 



We now turn our attention to convex sets for which the associated projection does not have a closed form, 
and we show that under some appropriate assumptions, it is possible to circumvent this difficulty. Let C 
denote such a nonempty closed convex subset of R M and assume that, for every y G R , 

yeC # h(y)=J2h {e) (y {l) )<V (23) 
t=i 

where r\ € K. Hereabove, the generic vector y has been decomposed into blocks of coordinates as follows 

?/=[(y (1) ) T ,.--,(y (i) ) T ] T (24) 

sizeM* 1 ) sizcMt-k) 

and, for every I G {1,...,L}, yW G R mW and is a function in T Q (R Mie) ) such that ri(dom/iW) ^ 0. 
Actually, any closed convex subset C can be expressed in this way by setting rj = 0, L = 1 and h = dc- 

By introducing now the auxiliary vector £ = (C^) 1<e<L G Condition ([25)1 can be equivalently 
rewritten at0 



;C W <r7, (25) 

,(Wg{i,...,l}) ^(y w )<C w - (26) 

Let us now introduce the closed half-space of R L defined as 

v = {<;eR L I ijc < (27) 

with !.£ = (1, . . . , 1) T G R L , and the closed convex set 

E={(y,C)eR M xR L | (We (y W ,C W ) € epi/i^}. (28) 

Then, Constraint ([23)1 means that ( € ^, whereas Constraint ([2^)1 is equivalent to (y, () G In other 
words, the constraint set C can be split into the two constraint sets V and E provided that L additional 
scalar variables (C )i<^<£ & r e introduced in the original optimization problem. Dealing with additional 
constraints in the original problem is not a problem for proximal splitting algorithms as far as the projections 
onto the associated constraint sets can be computed. An example of use of such an algorithm will be 
described in detail in Section @] . 

In the present case, the projection onto V is given by Proposition 12. 1[ whereas the projection onto E is 
given by 

(V(y, C) € K M x R L ) P E (y, C) = (p, 9) (29) 

where 9 — (0^)i<e<L, vector p G M M is blockwise decomposed as p = [(pW) T , . . . , (p^ L ') T ] T like in (|24l) . 
and 

(W G {1, . . . , L}) (pW 0«) = P epi/lW (yW, C W ). (30) 

So, the problem reduces to the lower-dimensional problem of the determination of the projection onto the 
convex subset epi using the results in Section 12.21 these 

projections can be shown to have a closed form expression in a number of cases of particular interest. Two 
examples are now given: 



2 Note that the inequality in i 125[ ) can be also replaced by an equality although it makes little difference in our approach. 
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• Euclidean norm functions of the form 



(Vy«GR MW ) ftW(yW) = rW||yW|| (31) 
where ^ G {1, . . . , L} and £ ]0, +oo[. 



Proposition 3.1. Assume that h^> is given by ([3T 
TTiera, for every (y w ,C W ) € M M<£> x R, 



•(0,0), */l|y (£) ll<-r«C w , 

^e Pi ^(y w ,C (£) ) = <|(yW,C w ), */lly w ll<S> ( 32 ) 

a w (y W ,r (£) ||y m ||), 



WfeereaW= l + (rW)^ ( 1 + 17^) 



The result of this proposition is a direct application of Proposition 12.41 and Corollary 12.71 As it will 
be shown in Sectional this result is useful to deal with multivariate sparsity constraints |74) or total 
variation bounds 75 , 76 , since such constraints typically involve a sum of functions like pip composed 
with linear operators corresponding to analysis transforms or gradient operators. 



Infinity norms defined as: for every £ <E {1, . . . , L} and = (y^' m - l )i< m <A/(«) € M M< , 

^(yW) = max{j^l|l<m<Mw|. 



(33) 



where {r^ m) ) 1<m<MW e]0,+c^[ 



mo 



Proposition 3.2. Assume that h™' is given by ([3"3"f . Lei (^' m 0i< m <MW ^ e a sequence of reals ob- 
tained by sorting (|y^' m ^|/T^ ,m )) 1 < m <j M -(«) in ascending order, and set i/^>°) — — oo and v^> m{ >+1 ) = 
+oo. Then, for every <W € R 7 (p^>W) = P opihW (yW C W ) *s suc/i ftaf p« = (p^^x^M^ , 

' y (Am) ) ^ |y(^,m)| < T {l,m)Q(i)^ 

T (««W, i/ y(^ ro ) > r^flW, (34) 



max(cW+E^^ m) (^) 2 ,0) 

^ " i + y M( Ur^Y ' (35) 



and m is the unique integer in {1, ... , M^> + 1} such that 
The proof of this proposition is given in Appendix [Fl 

Note that when r^ ,m ' = 1, function ftW in ([33]) reduces to the standard infinity norm || • H^. This 
proposition can be thus employed to efficient deal with £i j00 regularization which has attracted much 
interest recently [771 l78l [79] . 



S 



4 Experimental Results 



In this section, we provide numerical examples to illustrate the usefulness of the proposed epigraphical 
projection method. The presented experiments focus on applications involving projections onto ^i iP -balls 
where p £ {2, +00}. 

As already mentioned in Section |2.1[ various algorithms can be used to solve non-smooth convex opti- 
mization problems and would potentially benefit from the proposed epigraphical projection technique. In 
this work, we will employ a primal-dual algorithm, namely the Monotone+Lipschitz Forward Backward 
Forward (M+LFBF) algorithm, which was recently proposed in [68]. This proximal algorithm is able to 
address a wide class of convex optimization problems without requiring any matrix inversion and it offers a 
good performance and robustness to numerical errors. Its convergence is guaranteed (under weak conditions) 
and its structure makes it suitable for implementation on highly parallel architectures. 



4.0.1 Degradation model 

Set H — 1 N . Denote by x = {x^) 1<n< jj £ the signal of interest, and by z £ I 1 *' an observation vector 

such that z = Ax + b. It is assumed that A £ H NxN is a linear operator, and b £ M is a realization 
of a zero-mean white Gaussian noise. The recovery of x from the degraded observations is performed by 
following a variational approach which aims at solving the following problem 

L 

minimize \\Ax-z\\ 2 s. t. ^\\n t B t Fx\\ p <r}, (37) 

where (fi,~p) S K 2 with fi < ~p, rj is a real positive constant, and F £ M. KxN is the linear operator associated 

with an analysis transform. Furthermore, for every £ £ {1, . . . , L},Bg £ R M< xK is a block- selection linear 
operator which selects a block of data from its input vector U For every £ £ {1, . . . , L}, Vlt denotes an 
x diagonal matrix of real positive weights. 

The term \\Ax — z\\ 2 is the data fidelity corresponding to the minus log-likelihood of x. The bounds \i 
and ~p allow us to take into account the value range of each component of x. 

The second constraint involved in Problem (|37p conveys a smoothness condition. It exploits the fact 
that natural signals usually exhibit a smooth spatial behaviour, except around some locations (e.g. object 
edges in images) where discontinuities arise. The proposed constraint is based on a weighted £i iP -norm, 
where p £ {2, +00}. The associated block-sparsity measure extends many of the smoothness terms used 
in the literature. It reduces to the weighted £i-norm criterion found in [80] when each block reduces to a 
singleton (i.e. L — K, and, for every £ £ {1, . . . , L}, ~ 1 and Biy = y^')0 It captures the £x t 2 criteria 
present in [75] 81, 82, 83] when p — 2. It matches the t\ j00 criterion proposed in [78] when p = +00. 

Let us define A = [Bjn u Bjn L ] T and 

L 

C 2 = {y£R M I J2h { \<v} (38) 
e=i 

3 This means that there exist distinct indices mi, . . . ,m M (t) in {1, . . . , K} such that, for every y = {y^)i<k<K S R^, 
4 In this case, the constraint in I I37II reduces to a sum of absolute values. 
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where M = M^> + ■ ■ -+M^ L ' and the same decomposition as in (|2~41 is performed. Then, it can be observed 
that Problem (f57|) is a particular case of Problem 11.11 where S = 2, R — 1, g\ = \\ ■ — z\\ , T\ = A, Hi = I, 
C\ = \p, Jt] , i?2 = Ai* 1 , and C2 is the above £i. p -ball. 



4.0.2 Algorithmic solution 

The main difficulty in solving Problem (|37[) stems from the second constraint. The point is that most of the 
applicable algorithms require to compute the projection onto Ci. Specific numerical methods 84, 851 1751155] 
have been developed for this purpose. The aim of this section is to propose an alternative method based on 
the splitting principle presented in Section [3] So doing, the resulting problem can be efficiently addressed 
by proximal algorithms. The two possible approaches are now detailed. 

• Epigraphical method - The principle of this method is to decompose C2 into the closed half-space 
defined by (|2"7|) and the closed convex set defined by (j2"5j) with = \\-\\ p . The advantage of this 
decomposition is that the projections onto V and E have closed forms. Indeed, Py is given by 
Proposition [24j while Pg is given by Proposition 13. II for p = 2 and Proposition 13.21 for p = +00. We 
are then able to solve Problem (l3"Tj) by means of the M+LFBF algorithm. The associated iterations 
are given in Algorithm Q] 



Algorithm 1 Epigraphical method: M+LFBF for solving Problem (|37p 

Initialization 

6 = 2|jA|| 2 +max{||AF||,l} 
- ^]0,^l 
For i = 0, 1,.. . 

7, £ [e, i^S] 

(£[<], £W) = (xH,C w ) -7 l (2A T (AxW -z) + F T A T vM , i/W) 

(pW >p M ) = ( P *<8M),iV(?M)) 

(^[t] 5 p[i] ) = [IT y [*] ) + 7i ( AF x H , C W ) 

(aW,aH) = -7 l P E (7" 1 «;H, T- -1 ^) 

(t?M , ) = ( a [<] , Q W ) + 7l (AF p H , ) 

(«[»+!], = («[»] +??[«] -?W +57H) 

(JW,C H ) = (p w ,p H ) -7i(2A T (ApW _ 2 ) +F T A T aW,aW) 

(a.[«+il, £[*+*]) = ( X M -$[<] + +C H ) 



• Direct method - For completeness, we also consider the projection onto C2 with the algorithm in [84] 
when p — 20 or the iterative algorithm in j78j when p — +00H Then, M+LFBF can still be used to 
solve Problem ([37]) . 

Note that, according to the general result in |68[ Theorem 4.2], the sequence (a;M) generated by Algo- 
rithm [1] is guaranteed to converge to a (global) minimizer of Problem ([3"T]) . 

5 Codc available at www.cs.ubc.ca/~mpf/spgll 
6 Code available at www.lsi.upc.edu/~aquattoni 
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4.0.3 Smoothness constraint 



We consider two smoothness constraints based on recent Non-Local TV measures [5TJ [57] . They constitute 
particular cases of the one considered in (|37j) when L = N. They are described for 2D data in the following. 



• t 2 



£2-Non-Local TV - This constraint has the form 



N ,1/2 



53 ( 2 u^-x^f) < v , (39) 

1=1 nejVfCW, 

where A/i is the neighbourhood support of £, and is the set of positions n £ {1, . . . , iV} \ {£} located 
into a Q x Q window centered at i, where Q £ N is odd. This constraint is a particular case of the 
one considered in (|37|) where K = (Q 2 — 1)N and F is a concatenation of discrete difference operators 
F gi , ga with (gi,q 2 ) e {-(Q-l)/2,...,(Q-l)/2} 2 \{(0,0)}. More precisely for every (qi,q 2 ), F quq2 
is a 2D filter with impulse response: for every (ni,n2) € Z 2 , 

{1, if nx = n 2 = 0, 
-1, ifni=?i andn 2 = g 2l (40) 
0, otherwise. 

In addition, for every £ 6 {1, . . . , AT}, M^- 1 < Q 2 — 1, i?^ selects the components of Fx corresponding 
to differences (x^> — x^) n ^ t , and the positive weights ((jJe,n)neJ^e are gathered in the diagonal matrix 

t-oo-Non-Local TV- We consider the following constraint 

lv 

Y^max^nlxW -x^\} <r,. (41) 



i=l 



We proceed similarly to the previous constraint, except that the ^oo-norm is now substituted for the 
£ 2 -norm. 

Note that the classical isotropic TV constraint (designated by .£ 2 -TV in the following) constitutes a particular 
case of the £ 2 -NLTV one, where each neighbourhood Me only contains the horizontal/ vertical neighbouring 
pixels (M W = 2) and the weights are n = 1. Similarly, the ^-TV constraint is a special case of the 
^-NLTV one. 



4.0.4 Weight estimation and neighbourhood choice 

To set the weights, we got inspired from the Non-Local Means approach originally described in 88 . Here, for 
every £ € {1, . . . , TV} and n € Me, the weight ue.n depends on the similarity between patches built around 
the pixels £ and n of the image. Since our degradation process involves some missing data, a two-step 
approach has been adopted. In the first step, the -£ 2 -TV approach is used in order to obtain an estimate x 
of the target image. This estimate is subsequently used in the second step to compute the weights through 
a self-similarity measure, yielding 

uj e , n = uj e exp (-<5~ 2 \\B e F e x - B n F n x\\ 2 ^j , (42) 
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where S £ K \ {0}, loi g ]0, +oo[, Bi (resp. B n ) selects a Q x Q patch centered at position £ (resp. n) and 
.Fi (resp. F n ) is a linear processing of the image depending on the position £ (resp. n). The constant uii is 
set so as to normalize the weights (i.e. X^ n eA/" UJl - n = -*-)• 

The measure in (|42|) generalizes the one proposed in [88], which corresponds to the case when (resp. 
i 7 ^) reduces to a Gaussian function with mean £ (resp. n). In the present work, we consider the foveated self- 
similarity measure recently introduced in [89], due to its better performance in denoising. This approach 
can be derived from (j42j) by setting Ft (resp. F n ) to a set of low-pass Gaussian filters whose variances 
increase as the spatial distance from the patch center £ (resp. n) grows. 

For every I 6 {1, . . . , N}, the neighbourhood Mi is built according to the procedure described in [90] . 
In practice, we limit the size of the neighbourhood, so that < M. 

4.0.5 Simulation results 

The purpose of this section is twofold. First, the quality of images reconstructed with our variational 
approach is evaluated for different choices of regularization constraints and comparisons are made with a 
state-of-the-art method. Secondly, the convergence speed of the proposed epigraphical technique is evaluated 
w.r.t. the direct methods. 

In the following experiments, if not specified otherwise, the degradation matrix A is a decimated convo- 
lution which consists of a 3 x 3 uniform blur followed by a decimation which randomly deletes 60% of the 
pixels (N — 0.4 x N). The standard deviation of the additive white Gaussian noise is equal to o — 10. Since 
we deal with natural images, the data range bounds are \x = and ~p = 255. For the smoothness constraint, 

we set Q = 11, Q = 5, S = 35 and M = 14. 

Extensive tests have been carried out on several standard images of different sizes. The SNR and 
SSIM |91] results obtained by using the various previously introduced TV-like constraints are collected in 
Table [1] In addition, a comparison is performed between our method and the Gradient Projection for Sparse 
Reconstruction (GPSR) method [BO], which also relies on a variational approach. The constraint bound for 
both methods was hand-tuned in order to achieve the best SNR values. The best results are highlighted 
in bold. A visual comparison is made in Figure [1] where two representative images are displayed. These 
results demonstrate the interest of considering non-local smoothness measures. Non-Local TV with £1.2- 
norm indeed proves to be the most effective constraint with gains in SNR and SSIM (up to 1.82 dB and 
0.042) with respect to ^2-TV, which in turn outperforms GPSR. The better performance of NLTV seems 
to be related to its ability to better preserve edges and thin structures present in images. In terms of 
computational time, GPSR is about twice faster than ^2-NLTV. Our codes were developed in MATLAB^, 
the operators F and F T being implemented in C using mex files. 

In order to complete the analysis, we report in Figure[2]SNR/SSIM comparisons between ^2-NLTV and 
^2-TV for different blur and noise configurations. These plots show that ^2-NLTV provides better results 
regardless of the degradation conditions. 

In the second part of the section, we focus on the convergence speed of epigraph and direct methods. 
All the results refer to the culicoidae image cropped at 256 x 256 (N = 256 2 ), since a similar behaviour was 
observed for other images. The stopping criterion is set to: 

lla^+i] _ < io- 4 ||a;M||. For the £i iP -ball 
projectors needed by the direct method, we used the software publicly available on-line [84] [78] . 

7 R2011b version on an Intel Xcon CPU at 2.80 GHz and 8 GB of RAM. 
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Tabic 1: SNRdB and SSIM results of our method and GPSR (noise parameters: blur — 3 x 3, a = 10, 
decimation = 60%) 



SNR (dB) - SSIM 


N 


£2 


-TV 


p 


-TV 


£ 2 -NLTV 


^oo-NLTV 


GPSR 


CULICOIDAE 


256 2 


20.80 - 


- 0.855 


20.25 


- 0.853 


22.62 - 


0.897 


22.38 


- 0.897 


17.03 - 


0.738 


Lena 


256 2 


23.18 - 


- 0.783 


22.77 


- 0.769 


24.18 - 


0.812 


24.14 


- 0.812 


20.26 - 


0.678 


Cameraman 


256 2 


20.06 - 


- 0.774 


19.68 


- 0.755 


20.71 - 


0.801 


20.17 


- 0.743 


17.92 - 


0.673 


Boat 


256 2 


20.25 - 


- 0.739 


19.74 


0.718 


21.13 - 


0.770 


20.77 


- 0.741 


18.06 - 


0.649 


House 


256 2 


25.47- 


- 0.823 


24.70 


- 0.808 


26.31 


0.836 


25.87 


- 0.823 


22.14 - 


0.734 


Man 


256 2 


19.24 


- 0.725 


18.96 


0.714 


19.66 - 


0.741 


19.51 


- 0.736 


17.11 


0.629 


Peppers 


512 2 


23.69 - 


- 0.801 


23.25 


- 0.786 


24.80 - 


0.829 


24.45 


- 0.813 


21.94 - 


0.709 


Barbara 


512 2 


16.74 - 


- 0.653 


16.64 


- 0.642 


17.02 


0.673 


16.99 


- 0.652 


15.97- 


0.562 


Hill 


512 2 


22.18 - 


- 0.723 


21.89 


- 0.715 


22.55 


0.735 


22.43 


- 0.733 


20.21 - 


0.637 


CULICOIDAE 


1024 2 


20.84 - 


- 0.855 


20.49 


- 0.812 


23.25 - 


0.885 


22.57 


- 0.810 


17.19 - 


0.725 



H 



(a) Culicoidae. (b) Degraded 



(c) Zoom. 




(d) GPSR, 

SNR: 17.03 dB 



Zoom. (h) GPSR, 

SNR: 20.26 dB. 





*m jam jam jam - 

HHHH 



(i) I2-TV, (j) £oo-TV, (k) £ 2 -NLTV, (1) ^oo-NLTV, (m) £ 2 -TV, (n) f^-TV ', (o) £ 2 -NLTV, (p) £oo-NLTV, 
SNR: 20.80 dB. SNR: 20.25 dB. SNR: 22.62 dBBNR: 22.38 dB. SNR: 23.18 dB.SNR: 22.77 dB. SNR: 24.18 dBSNR: 24.14 dB. 

Figure 1: Image restoration examples (noise parameters: blur = 3 x 3, a = 10, decimation = 60%) 
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0.9 




0.8 




0.8 




0.7 




0.7 





15 



20 



25 



0.8 



0.9 



0.7 



0.9 



(a) SNR comparison (3x3 (b) SNR comparison (5x5 (c) SSIM comparison (3 X 3 (d) SSIM comparison (5x5 
blur). blur). blur). blur). 

Figure 2: SNR and SSIM values for £ 2 -NLTV (vertical axes) and £ 2 -TV (horizontal axes), for the culicoidae 
image. The plots show the results obtained for a € {5, 10, . . . , 50}, where lower SNR or SSIM values 
correspond to higher a values. No decimation is applied in this experiment. 
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• Total Variation - Tables [2] and [3] report a comparison between the direct and epigraphical methods for 
different values of r\ and for ^2-TV and ^oo-TV, respectively. For more readability, these are expressed 
as a multiplicative factor of the £ p -TV-semi-norm of the original image. It can be noticed that the 
parameter r\ influences both the quality of the results and the convergence speed. For the fi^-norm, 
the epigraphical method is 3.5 faster than the direct one. For the fi i00 -nonn, it is 65 times faster. 

Figs. 13-al and 13-bl show the relative error ||a;M — x' 00 ' H/H^ 00 ! || as a function of the computational time, 
where a^°°l denotes the solution computed after a large number of iterations (typically, 5000 iterations). 
The dashed line presents the results for the direct method while the solid line refers to the epigraphical 
one. These plots show that the epigraphical approach is faster despite it requires more iterations in 
order to converge. This can be explained by the computational cost of the sub-iterations required by 
the direct projections onto the ^i jP -ball. 



Table 2: Results for the ti-TN constraint and different values of r\ 





SNR (dB) - SSIM 


DIRECT 


EPIGRAPH. 




n 


SPEED UP 






# iter. sec. 


# iter. sec. 




0.6 


19.70 - 0.839 


116 5.74 


146 2.23 


2.58 


0.7 


20.62 - 0.862 


132 7.14 


151 2.29 


3.11 


0.8 


20.80 0.855 


160 9.17 


171 2.59 


3.54 


0.9 


20.45 - 0.826 


195 11.95 


196 2.98 


4.01 


Table 3: Results for the . 


(!oo-TV constraint and different values of i] 




SNR (dB) - SSIM 


DIRECT 


EPIGRAPH. 




n 


SPEED UP 






# iter. sec. 


# iter. sec. 




0.6 


19.58 - 0.839 


195 373.84 


230 6.21 


60.22 


0.7 


20.25 0.853 


206 417.71 


231 6.39 


65.34 


0.8 


20.24 - 0.835 


233 486.17 


245 6.86 


70.82 


0.9 


19.80 - 0.800 


263 551.03 


267 7.44 


74.06 



• £-2-Non-Local Total Variation - Table 2] reports the results for ^2-NLTV. Different combinations of 
neighbourhood size Q and bound value r\ are considered. To set the weights, the first TV estimate 
is computed with r\ = 0.8. The best SNR improvement (1.82 dB over ^2-TV) is observed when a 
relatively small neighbourhood is used (Q = 11). This behaviour was already observed in [90]. It may 
be related to the bias- variance trade-off commonly encountered in Non-Local filters [92j [93] . 

In Figure [3^cl a plot similar to those in Figs. 13-a l and !3-bl show the convergence profile. The epigraphical 
method requires about the same number of iterations as the direct one in order to converge. This 
results in a time reduction (1.8 times), as a single iteration of the epigraphical method is faster than 
one iteration of the direct method. 

• £oo-Non-Local Total Variation - Tableland Figure l3^dl show the results obtained with the £oo-NLTV 
constraint. Note that €oo-NLTV requires more iterations than .62-NLTV in order to converge. In what 
concerns the convergence speed, the epigraphical method is 19 times faster. 
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Table 4: Results for the ^2-NLTV constraint and some values of 77 and Q 



„. m OOT „ DIRECT EPIGRAPH. 
7] SNR (dB) - SSIM SPEED UP 







$ iter. sec. 


$ iter. 


sec. 








Neighbourhood sii^P" 


— 3 






0.7 


22.26 - 


0.895 70 5.82 


77 


3.07 


1.90 


0.8 


22.41 - 


0.893 72 6.39 


75 


3.00 


2.13 


0.9 


22.06 - 


0.875 88 7.97 


89 


3.58 


2.23 






Neighbourhood size: 


= 5 






0.7 


22.44 - 


0.900 70 7.13 


74 


4.37 


1.63 


0.8 


22.58 - 


0.898 72 7.85 


75 


4.44 


1.77 


0.9 


22.25 - 


0.880 87 9.71 


88 


5.24 


1.85 






Neighbourhood size: 


= 11 






0.7 


22.50 - 


0.901 72 7.52 


76 


4.51 


1.67 


0.8 


22.61 


0.897 75 8.21 


78 


4.64 


1.77 


0.9 


22.27 - 


0.879 89 9.83 


91 


5.41 


1.82 



Table 5: Results for the Iqq-NLTV constraint and some values of r t and g 

OOT ,. DIRECT EPIGRAPH. 

rj SNR (dB) SSIM speed UP 







# iter. sec. 


# iter. 


sec. 








Neighbourhood size: 


= 5 






0.6 


21.78 


- 0.889 229 766.21 


251 


25.61 


29.92 


0.7 


22.28 


- 0.894 213 653.31 


226 


22.87 


28.56 


0.8 


22.10 


- 0.876 207 596.90 


216 


21.89 


27.26 






Neighbourhood size: 


= 5 






0.6 


21.92 


- 0.894 231 914.03 


252 


43.32 


21.10 


0.7 


22.38 


0.897 219 749.32 


232 


39.73 


18.86 


0.8 


22.21 


- 0.880 212 673.60 


224 


38.05 


17.70 






Neighbourhood size: 


= 11 






0.6 


21.95 


- 0.894 236 940.38 


256 


43.75 


21.49 


0.7 


22.38 


- 0.897 222 760.49 


234 


39.84 


19.09 


0.8 


22.20 


- 0.879 216 689.07 


227 


38.80 


17.76 



5 Conclusions 



We have proposed a new epigraphical technique to solve constrained convex optimization problems with 
the help of proximal algorithms. In this paper (Part I), our attention has been turned to block-sparsity 
constraints based on weighted ^i^-norms with p E {2, +00}. The obtained results demonstrate the better 
performance of non-local measures in terms of image quality. Our results also show that the ^i^-norm 
has to be preferred over the ^i j00 -norm for image recovery problems. However, it would be interesting 
to consider alternative applications of .£i i00 -norms such as regression problems [9HI95]. Furthermore, the 
experimental part indicates that the epigraphical method converges faster than the approach based on 
the direct computation of the projections via standard iterative solutions. Full implementations in C of 
the proposed algorithms and parallelization of our codes should even allow us to accelerate them |96) . 
Note that, although the considered application involves two constraint sets, the proposed approach can 
handle an arbitrary number of convex constraints. In a companion paper (Part II), we will show that the 
epigraphical approach can also be used to develop approximation methods for dealing with more general 
convex constraints. 
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Figure 3: Comparison between epigraphical method (solid line) and direct method (dashed line): 
in dB vs time 
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A Proof of the Proposition I2.4I 



For every (y, () G % x R, let (p, 6>) = P op i V (y, C)- If ^0/) < C then p = y and = C = max{^(p),C}. In 
addition, 



(Vw G H) = y|| 2 + i(maxMp) - C,0}) 2 



1„ 



i 



max{v3(u) - C,0})' 



(43) 



which shows that (fT2")) holds. Let us now consider the case when <p(y) > C- From the definition of the 
projection, we get 

(p,6)= argmin \\u - y\\ 2 + (£ - () 2 . (44) 

(«,£)Gepi i/3 

From the Karush-Kuhn- Tucker theorem j97l Theorem 5.2] there exists a € [0, +oo[ such that 



(p,0)= argmin - y|| 2 + -(£ - C) 2 + a(<p(«) - 



where the Lagrange multiplier a is such that 



a(<p(p) - 0) = 0. 



(45) 



(46) 



Since the value a = is not allowable (since it would lead to p — y and 8 — £) , it can be deduced from the 
above equality that <p(p) = 9. In addition, differentiating the Lagrange functional in (|45[) w.r.t. £ yields 



tpip) = 9 = C + a>(. 



(47) 



8 By considering uq £ domip and £o < ¥>(wo)> the required qualification condition is obviously satisfied. 
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Hence, (p, 9) given by (|4"4l is such that 

p = argmin||u - y\\ 2 + (<p(u) ~ C) 2 (48) 
Mew 
v(u)>C 

9 = <p(p) = max{<p(j>), C}. (49) 

Furthermore, as <p(y) > £ we have 

inf ||u - yf = ||P lcv<c M ~y\\ 2 = mf ||u - y\\ 2 (50) 

¥>(«)<C »>(u)=C 

where we have used the fact that P\ C v <t - v (y) belongs to the boundary of lev<^ tp which is equal to 
{ii£M J <p(u) = since <p is lower-semicontinuous and dom^ is open |73[ Corollary 8.38]. We have 
then 

inf \\u- y\\ 2 = inf \\u - y\\ 2 
ip(u)<C ¥>(«)=£ 

> inf \\u-y\\ 2 + ^(u)-() 2 . (51) 

¥>(«)>£ 

Altogether, (gSJ and (J3TJ) lead to 



11 2 

p = argmin - y\\ 2 + -(max{^(w) - £ 0}) (52) 
ugh * *> 



which is equivalent to (|T2"j) since |( max{</? — £ 0}) 2 S T (^) 



B Proof of the Proposition 12.51 

Since (maxjip - £0}) 2 is an even function, proxi (max^-^o}) 2 1S an °dd function [331 Remark 4.1(h)]. In 
the following, we thus focus on the case when y £ ]0, +00 [. 

If C €] - 00, 0], then (max{(^ - £ 0}) 2 = (<p - () 2 . When = 1, (max{</> - C, 0}) 2 = r 2 (-) 2 - 2rC| • | + C 2 
and, from Example 12.31 it can be deduced that 

Proxitmax^-c.o}) 2 ^) =P r °x__ii 2 .|.| (j^j) = ^ ^ t2 max{y + r£0}. (53) 

When /3 > 1, (ip — Q 2 is differentiable and, according to ([8]). p = proxi^ max r ¥ ,_^ }) 2 (?/) i s uniquely defined 
as 

p-y + /3t P i3 - 1 (t P p -C) = (54) 
where, according to [SSI Corollary 2.5], p > 0. This allows us to deduce that p = \ . 

Let us now focus on the case when C e]0, +oo[. If y e]0, (C/t) 1//3 [, it can be deduced from Corollary 2.5], 
that p = proxi (roax{v _ f)0})2 (y) £ [0, (C/r) 1 /^. Since (V« G [0, (C/r) 1 /^) max{^(«) - C,0} = 0, © yields 
p = y. On the other hand if y > (C/t) 1 /^, as the proximity operator of a function from R to R is continuous 
and increasing O Proposition 2.4], p = proxi (max{v _ C0})2 (y) > proxi (max{l/ ,_ C0 }) 2 ((C/t) 1/ ' 3 ) = (C/t) 1//? - 
Since (maxjtp — ^, 0}) 2 is differentiable in this case, and (Vv > (C/t) 1 ^) (ma,x{(p(v) — £0}) 2 = (tv@ — () 2 , 
§8§ allows us to deduce that p is the unique value in [(C/ T ) 1 ^' 9 , +oo[ satisfying (|54"|) . It can be concluded 
that, when ( e]0, +oo[, (15) holds. 
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C Proof of the Proposition 12.61 



Let us notice that i(max{rd^, — (,0}) 2 = -0 ° dc where ip = i(max{r| • I' 3 — (, 0}) 2 . According to [351 
Proposition 2.7], for every y G %, 

(y, HyeC, 

P rox 0od c (y) = S p c(y), if y ^ C 1 and dc(y) < max 9-0(0), (55) 

[ay + (1 - a)P c (y), if G?c(y) > max #0(0) 

where a — - — *. y°^) . i n addition, we have 

^ ( 0) = { [TC '- rC] ' iK<0an d/3 = l, 
l{0}, otherwise, 

and, according to Proposition 12. 5[ when £ < 0, /3 = 1 and dc(y) < — T C P rox v> (^c(y)) = 0- These show 
that ([55]) reduces to (fP7]). 



D Proof of the Corollary E7T] 

As we have ip = rdc where C = {z}, the result follows from Proposition 12.61 and the expression of 
proxi( max { T |.|_£ }-)2 in Proposition 12.51 



E Proof of the Proposition 12.81 

The function ip belongs to To(R) since for every m G {1,...,M}, v i-> max{r' m '(i/' m ' — v),0} is hnitc 
convex and (-) 2 is finite convex and increasing on [0,+oo[. In addition, ip is differentiable and it is such 
that, for every v G R and every k G {1, . . . , M + 1}, 

m—l m=k 

For every y G K, as p — prox (j/) is characterized by (|HJ|, there exists to G {1, ...,M + 1} such that 
yCm-l) < p < v {m) and 

y-P= £(ri m) ) 2 (p^ (m) ) + E_( T i m) ) 2 ^-^ (m) )- ( 58 ) 

m— 1 m— m 

This yields and we have : <4> i/™ -1 ) < p < The uniqueness of to G {1, . . . , M + 1} satisfying 

this inequality follows from the uniqueness of prox (y). 
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F Proof of the Proposition 13.21 



For every (yW,£M) e R M x R, in order to determine P cpihW (y (£) , C W ) 

we have to find 

min f(^)-C^) 2 + min || p M_ y WfY (59 ) 

|p(',M W )|< T (<.J< W ) g<«) 

For every #w G [0, +oo[, the inner minimization is achieved when, for every j G {1, . . . , M^'}, p^- m "> is the 
projection of y(*> m ) onto [— r^' m ) $W : r^' m > d^], which is given by (j3"4")l . Then, the problem reduces to 

minimize ((0<'> - C W ) 2 + V (max{|y(^™)| - t^9^,0}) 2 + l [0 , +oo[ (6^)) (60) 

m— 1 

which is also equivalent to calculate prox v+t[o +oo[ (C )i where ip is such that 

M< f > 

(W G R) = 2 H (max{r ( ^ m) (i/^' ro ^ - w), 0}) 2 . (61) 

m— 1 

By using now the fact that prox^ + ( = -P[o.+oo[°P rox tp ([Ml Proposition 12]) and by invoking Proposition 
12.81 the expression of the optimal solution in (|35]) follows. 
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